#pragma once

#include "GhostFiller.h"

/**
 * @file FVOperator.h
 * @author xingyifan
 * @date 2022-08-31 22:43
 *
 * @description: 网格插值和约束过程
 */

template <int Dim>
class FVGridOperator
{
public:
    FVGridOperator()
    {
        // 计算权的分布
        Tensor<1> baseWeight({1, 2, 1});
        weight_ = TensorOp::power<Dim>(baseWeight);
    }

    /**
     * @file TensorOperator.h
     * @author xingyifan
     * @date 2022-07-26 00:19
     *
     * @description: 满权算子
     */

    Grid<Dim> FullWeighting(Grid<Dim> &t)
    {
        // 粗网格规模
        int n2 = (t.width() - 1) / 2;

        // 粗网格
        Grid<Dim> tj(n2 + 1);

        do
        {
            // 将粗网格的指标发送到细网格上的对应位置
            tj.fine(t);

            // 不在边界上才进行加权，否则直接赋值
            if (tj.findBound() == -1)
            {
                // 获取周围一圈的值
                Tensor<Dim> value = t.cell();

                // 计算加权平均
                tj.value() = TensorOp::dot(weight_, value) / pow(4, Dim);
            }
            else
            {
                tj.value() = t.value();
            }
        } while (tj++);

        // 刷新指标
        t.fresh();

        return tj;
    }

    /**
     * @file FVOperator.h
     * @author xingyifan
     * @date 2022-08-31 22:44
     *
     * @description: 线性插值
     */

    Grid<Dim> Linear(Grid<Dim> &t)
    {
        // 细网格规模
        int n2 = (t.width() - 1) * 2;

        // 细网格
        Grid<Dim> tj(n2 + 1, 0);

        // 先插不动的点
        do
        {
            // 获得细网格上对应的指标
            t.fine(tj);
            tj.value() = t.value();
        } while (t++);

        // 刷新指标
        tj.fresh();

        // 按照维数插值
        for (int d = 0; d < Dim; d++)
        {
            // 沿此方向的步长
            int step = IntOp::power(n2 + 1, d);

            do
            {
                // 从 d 开始，如果该方向上的下标不是偶数，就增加 1
                for (int d1 = d; d1 < Dim; d1++)
                {
                    if (tj.index(d1) % 2 == 1)
                    {
                        tj.index(d1) += 1;
                    }
                }
                // 当沿着此方向不是边界点，就对左右两个点平均插值
                if (!tj.isBound(d))
                {
                    tj.value(step) = (tj.value() + tj.value(2 * step)) / 2;
                    tj.value(-step) = (tj.value() + tj.value(-2 * step)) / 2;
                }
            } while (tj++);
        }

        return tj;
    }

private:
    Tensor<Dim> weight_; // 权
};

/**
 * @file FVOperator.h
 * @author xingyifan
 * @date 2022-08-29 14:49
 *
 * @description: Laplacian
 */

template <int Dim>
class Laplacian
{
public:
    /**
     * @file FVOperator.h
     * @author xingyifan
     * @date 2022-08-29 21:29
     *
     * @description: 作为 Av=f 中 A 的形式，需要在 f 中补足边界条件
     */

    Tensor<Dim> operator()(Tensor<Dim> &v, GhostFiller<Dim> &ghostfiller, Real t = -1)
    {
        // 网格规模
        int n = v.width();

        // 存放结果
        Tensor<Dim> res(n);

        // 每段长度
        Real h = 1.0 / n;

        // 中心系数，它是 A[i,i] 元素的值
        Real center = -30 * Dim / (12 * h * h);

        // 遍历每一个点
        do
        {
            // 填充 ghost cell
            Vector ghost = ghostfiller(v);

            // 如果 t 不为 -1 就增加边界条件
            if (t >= 0)
            {
                ghost += ghostfiller.bdCond(v, t);
            }

            // 记录总和
            Real sum = 0;

            // 加权迭代
            for (int d = 0; d < Dim; d++)
            {
                sum += (-ghost[5 * d] + 16 * ghost[5 * d + 1] +
                        16 * ghost[5 * d + 3] - ghost[5 * d + 4]) /
                       (12 * h * h);
            }

            res.value() = sum + v.value() * center;
        } while (v++, res++);

        return res;
    }

    /**
     * @file FVOperator.h
     * @author xingyifan
     * @date 2022-08-29 21:31
     *
     * @description: 为 f 补充边界条件
     */

    void supplement(Tensor<Dim> &f, GhostFiller<Dim> &ghostfiller, Real t, Real coeff)
    {
        // 网格规模
        int n = f.width();

        // 每段长度
        Real h = 1.0 / n;

        // 遍历每一个点
        do
        {
            // 获得边界条件
            Vector bdvalue = ghostfiller.bdCond(f, t);

            // 记录总和
            Real sum = 0;

            // 加权迭代
            for (int d = 0; d < Dim; d++)
            {
                sum += (-bdvalue[5 * d] + 16 * bdvalue[5 * d + 1] +
                        16 * bdvalue[5 * d + 3] - bdvalue[5 * d + 4]) /
                       (12 * h * h);
            }

            // 修正值
            f.value() += sum * coeff;
        } while (f++);
    }
};

/**
 * @file FVOperator.h
 * @author xingyifan
 * @date 2022-08-31 17:29
 *
 * @description: Gradient
 */

template <int Dim>
class Gradient
{
public:
    GridX<Dim> operator()(Tensor<Dim> &phi, GhostFiller<Dim> &ghostfiller)
    {
        // 系统规模
        int n = phi.width();

        // 网格间距
        Real h = 1.0 / n;

        // 存放结果
        GridX<Dim> G(Dim, Grid<Dim>(n));

        // 产生一个单元区域
        Domain<Dim> domain;

        // 遍历每一个单元
        do
        {
            // 在每个方向上计算梯度，然后相加
            for (int d = 0; d < Dim; d++)
            {
                // 此方向上的填充
                Vector cell(5);

                // 获得 ghost cell
                ghostfiller(phi, cell.ptr(), d);

                // 此方向上的梯度
                G[d].value() = (-cell[4] + 8 * cell[3] - 8 * cell[1] + cell[0]) / (12 * h);

                // 速度指标移动
                G[d]++;
            }
        } while (phi++);

        return G;
    }
};

/**
 * @file FVOperator.h
 * @author xingyifan
 * @date 2022-08-29 20:43
 *
 * @description: Divergence
 */

template <int Dim>
class Divergence
{
public:
    /**
     * @file FVOperator.h
     * @author xingyifan
     * @date 2022-08-31 17:39
     *
     * @description: 标准散度填充
     */

    Tensor<Dim> operator()(GridX<Dim> &u, GhostFiller<Dim> &ghostfiller)
    {
        // 网格规模
        int n = u.front().width();

        // 网格间距
        Real h = 1.0 / n;

        // 存放结果
        Tensor<Dim> Du(n);

        // 产生一个单元区域
        Domain<Dim> domain;

        // 遍历每一个单元
        do
        {
            // 散度
            Real D = 0;

            // 在每个方向上计算梯度，然后相加
            for (int d = 0; d < Dim; d++)
            {
                // 此方向上的速度填充
                Vector cell(5);

                // 获得 ghost cell
                ghostfiller(u[d], cell.ptr(), d);

                // 此方向上的梯度
                D += (-cell[4] + 8 * cell[3] - 8 * cell[1] + cell[0]) / (12 * h);

                // 速度指标移动
                u[d]++;
            }

            Du.value() = D;
        } while (Du++);

        return Du;
    }

    /**
     * @file FVOperator.h
     * @author xingyifan
     * @date 2022-09-01 06:14
     *
     * @description: 计算 <ud phi>
     */

    Tensor<Dim> operator()(Tensor<Dim> &phi, GridX<Dim> &u, GhostFiller<Dim> &ghostfiller, Real t)
    {
        // 网格规模
        int n = phi.width();

        // 网格间距
        Real h = 1.0 / n;

        // 存放结果
        Tensor<Dim> Ladvf(n);

        // 产生一个单元区域
        Domain<Dim> domain;

        // 遍历每一个单元
        do
        {
            // 散度
            Real D = 0;

            // 在每个方向上计算梯度，然后相加
            for (int d = 0; d < Dim; d++)
            {
                // 如果是周期边界
                if (ghostfiller.type_[2 * d] == 2)
                {
                    D += (product(u[d], phi, ghostfiller, d, 1) - product(u[d], phi, ghostfiller, d, -1)) / h;
                }
                // 不是周期边界，假定 u 是常数
                else
                {
                    // 此方向上的 ghost cell
                    Vector phicell(5);

                    // 填充两个 ghost
                    ghostfiller(phi, phicell.ptr(), d);

                    // 计算边界
                    Vector phibd(5, 0);
                    ghostfiller.bdCond(phi, phibd.ptr(), d, t);

                    // 合并两个 ghost cell
                    phicell += phibd;

                    // 此方向上的梯度
                    D += u[d].value() * (-phicell[4] + 8 * phicell[3] - 8 * phicell[1] + phicell[0]) / (12 * h);
                }

                // 速度指标移动
                u[d]++;
            }

            Ladvf.value() = D;
        } while (phi++, Ladvf++);

        return Ladvf;
    }

    /**
     * @file FVOperator.h
     * @author xingyifan
     * @date 2022-08-31 22:46
     *
     * @description: 直接返回 D<uu>
     */

    GridX<Dim> operator()(GridX<Dim> &u)
    {
        // 网格宽度
        int n = u.front().width();

        // 宽度
        Real h = 1.0 / n;

        // 存放结果
        GridX<Dim> Duu(Dim, Grid<Dim>(n));

        // 创建填充器
        GhostFiller<Dim> ghostfiller;

        // 退出检测
        bool check = true;

        // 遍历每一个单元
        do
        {
            // 存放结果
            Vector sum(Dim, 0);

            // 遍历每个方向
            for (int d = 0; d < Dim; d++)
            {
                // 遍历每个方向
                for (int d1 = 0; d1 < Dim; d1++)
                {
                    sum[d] += (product(u[d1], u[d], ghostfiller, d1, 1) -
                               product(u[d1], u[d], ghostfiller, d1, -1)) /
                              h;
                }
                Duu[d].value() = sum[d];
            }

            // 每层的指标都递增
            for (int d = 0; d < Dim; d++)
            {
                // 当到达边界就退出循环
                check = u[d]++;
                Duu[d]++;
            }
        } while (check);

        return Duu;
    }

private:
    /**
     * @file FVOperator.h
     * @author xingyifan
     * @date 2022-08-31 11:31
     *
     * @description: F<phi,psi>_0.5ed = <phi>_0.5ed <psi>_0.5ed + h^2/12 sum_{d'\neq d}(Gd'\phi)_0.5ed(Gd'\psi)_0.5ed
     */

    Real product(Grid<Dim> &phi, Grid<Dim> &psi, GhostFiller<Dim> &ghostfiller, int d, int sign)
    {
        // 系统规模
        int n = phi.width();

        // 宽度
        Real h = 1.0 / n;

        // 计算 <phi>_0.5ed <psi>_0.5ed
        Real phiface = faceAverage(phi, ghostfiller, d, sign);
        Real psiface = faceAverage(psi, ghostfiller, d, sign);

        Real Gd = 0;
        for (int d1 = 0; d1 < Dim; d1++)
        {
            // 计算 G<d1 phi>i+0.5ed
            if (d1 != d)
            {
                Gd += gradient(phi, ghostfiller, d, d1, sign) * gradient(psi, ghostfiller, d, d1, sign) * h * h / 12;
            }
        }

        return phiface * psiface + Gd;
    }

    /**
     * @file FVOperator.h
     * @author xingyifan
     * @date 2022-08-31 13:44
     *
     * @description: discrete gradient operator in the transverse directions
     */

    Real gradient(Grid<Dim> &phi, GhostFiller<Dim> &ghostfiller, int d, int d1, int sign)
    {
        // 系统规模
        int n = phi.width();

        // 宽度
        Real h = 1.0 / n;

        // 记录 d1 方向的指标
        int tind = phi.index(d1);

        // 向 d1 正向移动
        phi.forward(d1);

        // 计算面积分
        Real phiface1 = faceAverage(phi, ghostfiller, d, sign);

        // 恢复指标
        phi.index(d1) = tind;

        // 向 d1 反向移动
        phi.backward(d1);

        // 计算面积分
        Real phiface2 = faceAverage(phi, ghostfiller, d, sign);

        // 恢复指标
        phi.index(d1) = tind;

        return (phiface1 - phiface2) / (2 * h);
    }

    /**
     * @file FVOperator.h
     * @author xingyifan
     * @date 2022-08-31 13:49
     *
     * @description: 计算面积分
     */

    Real faceAverage(Grid<Dim> &phi, GhostFiller<Dim> &ghostfiller, int d, int sign)
    {
        // 填充 ghost
        Vector phii(5);

        // 获得 d 维上的 ghost
        ghostfiller(phi, phii.ptr(), d);

        // 填充中间值
        phii[2] = phi.value();

        // 计算 <phi>_0.5ed <psi>_0.5ed
        return (phii[2] + phii[2 + sign]) * (7.0 / 12) - (phii[2 - sign] + phii[2 + 2 * sign]) * (1.0 / 12);
    }
};